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Abstract 

A standard way to solve linear algebraic systems Au = f, (*) with ill-conditioned 
matrices A is to use variational regularization. This leads to solving the equation 
(A*A + aI)u = A*fs, where a is a regularization parameter, and f$ are noisy data, 
11/ — fs\\ < S. Numerically it requires to calculate products of matrices A* A and 
inversion of the matrix A* A + al which is also ill-conditioned if a > is small. We 
propose a new method for solving (*) stably, given noisy data f$. This method, the 
DSM (Dynamical Systems Method) is developed in this paper for selfadjoint A. It 
consists in solving a Cauchy problem for systems of ordinary differential equations. 



1 Introduction 

Consider a linear algebraic system 

A*Au = f, fER(A), (1) 

where A is a linear operator in n— dimensional Euclidean space, A = A*, R(A) is the 
range of A, and N := {u : Au = 0} is the null-space of A. Let k(A) : = | \A\ 1 1 {A' 1 ] | 
denote the condition number of A. If A is singular, i.e., N is not trivial, then we set 
k(A) = oo. Problem (1.1) is called ill-conditioned if k(A) >> 1. In this case small 
perturbations of / may lead to large perturbations of the solution u, so problem (1.1) 
is ill-posed. Such problems are often solved by variational regularization. This method 
consists in finding global minimizer of the functional F(u) = \ \Au — fs\\ 2 + a\\u\\ 2 , where 
a = const > is a regularization parameter, and f$ are noisy data, ||/ — fg\\ < 5. 
The global minimizer of the quadratic functional F is the unique solution to the linear 
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algebraic system (A*Au + al)u a j = A* f$, where / is the unit matrix. This system 
has a unique solution u a ^ = (A*Au + aI)~ 1 A*f$. Calculation of the matrix A* A requires 
multiplication of two matrices. For large n this is a time-consuming operation. Condition 
number of the matrix A* A is k 2 (A), so it is much larger than k(A) since k(A) >> 1. If 
o is small, then the condition number of the matrix A* A + al is also large. Therefore, 
inversion of the matrix A* A + al is numerically difficult if a is small. An additional 
difficulty consists in choosing the regularization parameter a as a function of 5 in such a 
way that the element us := « a (<5),<5 would converge to a solution of (1.1) as 5 — > 0. There 
are a priori and a posteriori methods for choosing such a(5). The theory of variational 
regularization is presented in many books and papers (see e.g., [1], Chapter 2). 

Our goal is to develop a new method for for solving ill-conditioned problems (1.1). 
This method, which we call the DSM (Dynamical Systems Method), does not require 
inversion of matrices and their multiplications. It requires solving a Cauchy problem 
for a system of ordinary diferential equations. The theoretical development of DSM is 
presented in [2], [3]. The author hopes that DSM will be an efficient numerical method 
for solving ill-conditioned linear algebraic systems. The assumption that A is selfadjoint, 
which is used for simplicity in this paper, can be relaxed. However, all generalizations 
will be considered elsewhere. Here we concentrate on the assumptions leading to the 
simplest arguments. 

The idea of our method is simple: it is based on the formula 

B-\e Bt -I) = f e Bs ds, (2) 
Jo 

where B is a linear boundedly invertible operator. If 

lim ||e Bt || = 0, (3) 

t—*oo 

then 

B' 1 = - lim / e Bs ds. (4) 

The integral J* e B ^~ s ^dsf is the solution to the Cauchy problem 

du 

u = Bu + f, u(0) = 0, u = —. (5) 

UjV 

Therefore, one can calculate the inverse of an operator satisfying condition (3) by solving 
a Cauchy problem. We want to calculate A' 1 /. Let us take B = i(A + ial), where / is 
the identity operator and a > is a parameter which we take to zero later. Condition 
(3) is satisfied if A = A* and a > because under these assumptions ||e Jj4 *|| = 1 and 
lim^oo e~ at = 0. We write A + ia in place of A + ial below. 
Consider the problem: 

u a = i(A + ia)u a + /, u a (0) = 0. (6) 
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Its unique solution is u a = J* * e l ( A + m )( t - s ) 'dsf \ Our results are the following two theorems. 
Theorem 1. One has 

—i lim lim u a (t) = y, (7) 

where y is the minimal-norm solution to (1.1). 
Note that if 

v a = i(A + ia)v a -if, v a (0) = 0, 

then lim a ^o hmt^oo v a (t) = y. 

If fs is given in place of /, || / — fs\ \ < 6, then one solves the problem: 

u a , s = i(A + ia)u a , s + fs, ii a ,*(0) = 0. (8) 
Theorem 2. If t$ and a = a(S) are such that 

limt 5 = oo, lima(5) = 0, lim— — = 0, lim a(8)tx — oo, (9) 
,5^0 <5^0 <5^0 a(S) s^o 



then 



lim \ \us — iy\ \ — 0, (10) 

8^0 



where u s := u a ( S ),s- 

In the next Section proofs are given. 



2 Proofs 

Proof of Theorem 1. By the argument in Section 1, we have 

Ua (t) = [- i (A + za)]- l (I-e« A+ ^ t )f, (11) 
and lim^oo ||e^ A+m ^|| = 0. Thus, lim^oo u a (t) = i(A + ia)~ l f . Consequently, 

-i\m\i(A + iaY 1 f = lim (A + ia)~ l Ay = y, 

as claimed in Theorem 1. Let us explain the last step. Using the spectral theorem for 
the selfadjoint operator A, one gets 

lim WiA+iaY 1 Ay -y\\ 2 = lim f |s(s+m) _1 -l\ 2 d(E s y,y) = lim f a 2 d{E s y,y) = 0. 

a->0 a^0 J R a^0 J R a z + S z 

(12) 

Here we have used the assumption y~LN, which implies that lim 0< ; ) ^o f_i,d(E s y,y) = 0. 
This relation allows one to pass to the limit a — > under the integral sign in (12). 
Theorem 1 is proved. □ 
Remark 1. In our case the operator A is bounded, so the integration in (12) is taken 
over a finite interval [— \\A\\, \ \A\\]. Moreover, in a finite-dimensional space the spectrum 
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of the operator A consists of finitely many eigenvalues, and the integral in (12) reduces to 
a finite sum. Our proof of Theorem 1 allows one to use it in infinite-dimensional spaces. 
Proof of Theorem 2. One has 

\\u s - 2/| | < \\u a{5) - 2/| | + || /"^ e« A+ia ^- s Us(fs - f)\\. (13) 

Jo 

The argument given in the proof of Theorem 1 shows that 

lim \\u a( s)(t s ) - 2/ 1 1 = (14) 

o— >0 



lim a(5) = 0, lim fx = oo, and lim a(5)ts = oo. 

5^0 S->0 6^0 



provided that 

One estimates the integral 

J HI f 6 e i{A+ia){t - s) ds(fs-f)\\< 
Jo 



(15) 



Thus 



\\u s -iy\\<- + o(l), (16) 

CL 

where the term o(l) comes from (14). Therefore, if lim^o — then lim^o I \us — y\ \ = 
0, and Theorem 2 is proved. □ 
Remark 2. The assumption A = A* allows one to give a short and simple proof of 
Theorems 1,2. However, this assumption can be replaced by more general assumptions. 
For example, one can assume that A has Jordan chains of length one only. In other words, 
that the resolvent of A has only simple poles. Under this more general assumption the 
Cauchy problem we have used should be also modified, in general. If A = A*, then our 
proofs of theorems 1 and 2 remain valid for operator A in a Hilbert space, and not only 
in a finite-dimensional space. 
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